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Abstract 

In the classical s-t network reliability problem a fixed network G is 
given including two designated vertices s and t (called terminals). The 
edges are subject to independent random failure, and the task is to com- 
pute the probability that s and t are connected in the resulting network, 
which is known to be #P-complete. In this paper we are interested in ap- 
proximating the s-t reliability in case of a directed acyclic original network 
G. We introduce and analyze an specialized version of a Monte-Carlo al- 
gorithm given by Karp and Luby. For the case of uniform edge failure 
probabilities, we give a worst-case bound on the number of samples that 
have to be drawn to obtain an e — 5 approximation, being sharper than 
the original upper bound. We also derive a variance reduction of the es- 
timator which reduces the expected number of iterations to perform to 
achieve the desired accuracy when applied in conjunction with difi'erent 
stopping rules. Initial computational results on two types of random net- 
works (directed acyclic Delaunay graphs and a slightly modified version 
of a classical random graph) with up to one million vertices are presented. 
These results show the advantage of the introduced Monte-Carlo approach 
compared to direct simulation when small reliabilities have to be estimated 
and demonstrate its applicability on large-scale instances. 

1 Introduction 

In the classical s-t network reliability problem, a fixed graph G with two special 
vertices s and t is given whose edges fail (disappear) independently of each other 
with some given probability. The task is to determine the probability that s and 
t are still connected in the resulting network after edge failures. A famous related 
problem is the all-terminal reliability problem, where the goal is to determine 
the probability that all vertices are still connected to each other after edge 
failures (for further information on various reliability problems see Both 
problems are known to be computationally hard (#P-complete) even for very 



restricted classes of graphs G [inillS]- In particular, the s-t reliability problem 
remains hard in the case when G is a directed acyclic planar graph of maximal 
degree three |12| . Therefore we are interested in finding approximations. 

Randomized algorithms respecting some relative error bound with high prob- 
ability have shown to be interesting approaches for different reliability problems 
|nj[71[5]. As we are interested in relative error bounds, there is a significant dif- 
ference between the estimation of the failure probability of the network and esti- 
mating the probability that the network is intact after edge failures. This results 
from the fact that estimating small probabilities by sampling is in general much 
more difficult than estimating large probabilities. Moreover, the techniques used 
for the estimation of reliability values often differ significantly from those used 
for the estimation of the failure probability. Most of the literature concentrates 
on estimating the failure probability because in typical applications, as com- 
munication or electrical networks, we have a highly reliable network and are 
interested in accurately estimating the probability of rare failures. 

In this paper, we concentrate on estimating the probability that the network 
is intact, i.e., the probability that there is a path from s to t after edge failures. 
This is motivated by certain models of spreading processes on networks, such as 
disease spreading, which can be mapped onto reliability problems |4l [T4 l [T7 l ITO l 
[TT| . In this context, a path from s to i represents the spread of a disease from 
s to t, and we are interested in estimating the probability of the rare event that 
the disease spreads over a long distance from s to t . 

The basis of our approach is a method presented by Karp and Luby |9j for 
estimating the failure probability in an s-t reliability problem when the graph G 
is planar, which can be easily adapted for estimating the probability of connect- 
edness of s and t on an arbitrary graph G. To be efficiently applicable, however, 
the underlying graph G needs to fulfill some additional properties, such as very 
low intactness probabilities of the edges and low vertex degrees. Furthermore, 
a computationally expensive preprocessing phase is required, which computes 
different quantities needed to efficiently sample from the proposed sample space. 
This makes the method not suited for application on large networks. We show 
in this paper that in the case of a directed acyclic network G, some of these 
problems can be circumvented, and the resulting algorithm can be applied to 
large-scale instances. The simplifications are due to the fact that generating an 
s-t path uniformly at random in an acyclic graph can be performed in linear 
time (see Section [5]) . Except for very restricted classes of initial networks such 
as series-parallel networks (see |15| for further information), no practically use- 
ful methods for estimating small s-t reliabilities on large networks are known. 
Considering directed acyclic graphs can thus be seen as a natural next step. 

The paper is organized as follows. We begin with some preliminaries in 
Section |2] In Section [3] the direct Monte-Carlo approach is described and its 
efficiency is analyzed. Section [4] then presents our algorithm for the estimation 
of small reliabilities in directed acyclic graphs up to some technical details that 
are explained in Section [sj In Section |6] we discuss how many samples have to 
be drawn in our algorithm to obtain a good estimate of the reliability with high 
probability. Section [7] contains computational results on two types of random 
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directed acyclic networks comparing the direct Monte-Carlo approach with our 
algorithm and demonstrating the applicability of our algorithm on large-scale 
instances. 



2 Preliminaries 

Let G — {V,E,p) be a directed acyclic network, where 
9 V is the set of vertices (with \V\ — n), 
9 E is the set of edges (with \E\ = m), and 

• p : E — > [0, 1] is a function associating a failure probability with every 
edge of the network. 

We call a network G satisfying these properties an acyclic reliability network. 
Furthermore, we fix two special vertices s,t ^V. Every edge e £ E fails with 
probability p{e) independently of the others. Let G" be the resulting graph over 
the intact edges after the realization of the edge failures. We say that G' is 
intact if it contains a path from s to t, otherwise we say that G' is in a failed 
state. Finally, the s-t reliability RELs.t{G) of the network G is defined as the 
probability of G' being intact. 

Let q : E — > [0, 1] be the function which associates with every edge its 
probability of being intact, i.e., q{e) = 1 — p{e) Ve € E. It is sometimes easier 
to look at our reliability model in a different way where edges appear rather 
than disappear. In this model we would begin with an empty network and flip 
a biased coin for all potential edges e G E to determine whether they appear. 

We are mainly interested in e — S approximations for RELg tiG), i.e., algo- 
rithms returning an estimate of RELs^t{G) accurate to within a relative error 
of e with probability at least \ ~ 5. To determine how many samples we have 
to draw in a Monte-Carlo algorithm to obtain an e — 5 approximation we often 
refer to the Generalized Zero-One Estimator Theorem introduced by Dagum et 
al. in [5]. The theorem is repeated below. 

Theorem 1 (Generalized Zero-One Estimator Theorem). Let X he a random 
variable taking values in [0, 1] and let Xi,X2, ■ ■ ■ ,Xn denote independent ran- 
dom variables distributed identically to X. If e < I and N > 4(e — 2)ln(2/(5) • 
max{Var[X],eE[X]}/(eE[X])2 th 



en 
1 ^ 

(l-e)E[X]<-5]X, <(l + e)E[X] 



>l-S . 



We typically use the theorem in the following form. If > 4(e — 2) ln{2/d) ■ 
(l/e^E[X]) then (X^i^i^O/-^ is an e — (5 approximation for E[X] (using the 
fact that a random variable X taking values in [0, 1] satisfies E[X] > Var[X]). 



3 A direct Monte- Carlo approach 

In this section we consider a simple Monte-Carlo approach and show that it is 
efficient for sufficiently large values of RELs^t{G), but inefficient for the estima- 
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tion of small reliabilities. In this approach we simply flip a biased coin for every 
edge e G E and observe whether s and t are connected in the resulting graph. 
Let X be the random variable corresponding to this approach where X = lif the 
resulting network is intact and otherwise. The random variable X has thus a 
Bernoulli distribution with parameter RELg tiG) and the reliability is estimated 
without bias by generating N independent realizations of X and returning their 
empirical mean which we denote by Y^r. A central question when using this 
method is how large N has to be chosen to obtain an e — 5 approximation. A 
direct application of Theorem [l] gives the following: 

Theorem 2. Y^r is an e — 5 approximation of RELs^t{G) if N satisfies 



When RELs^t{G) is bounded below by l/poly{n,m), is an FPRAS (e — 
6 approximation algorithm with running time bounded by a polynomial in 
1/e, log(l /S) and the input size) . The difficult case is the estimation of small val- 
ues, i.e., small reliabilities RELg^tiG). This problem motivated the construction 
of the algorithm to be presented next. 

4 Monte-Carlo method for estimating small re- 
liabilities 

The backbone of our algorithm is an adaption of the Monte-Carlo method pre- 
sented in [S]. Our algorithm exploits that, in a directed acyclic network, we can 
easily (in linear time) calculate the mean number of intact paths from s to t after 
the edge failures. This value is normally a good estimate for the reliability in 
highly unreliable networks as in such networks an intact state typically contains 
only a few paths from s to t. Using a Monte-Carlo approach, our algorithm then 
estimates the ratio between RELs t{G) and the mean number of intact paths 
from s to i after edge failures. Multiplying this estimate with the mean number 
of intact paths from s to t yields finally an estimate for the s-t reliability of the 
network. 

Note that the ratio between the mean number of intact paths from s to t after 
edges failures and the reliability to estimate, i.e., the reciprocal of the value we 
estimate in our algorithm, is exactly the mean number of intact paths from s to t 
after edge failures conditioned on the event that the network will be intact after 
the failure process. One of the main problems for the development of methods 
for directly estimating this ratio is the difficulty of choosing a sample out of 
the pool of intact states with probability proportional to its real appearance 
probability. In fact such a sampling procedure could easily be transformed 
into an FPRAS for the estimation of RELs,t{G) using techniques presented by 
Jerrum et al. 
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4.1 Notation 



Let H ^ {h C E} be the set of all possible states of the network G after the 
realization of the edge failures, where a state h G H represents the collection 
of intact edges. Furthermore, let A C H be the set of all intact states. For 
every state h £ H we denote by ■w{h) the probability that h occurs after the 
edge failure process, i.e., w{h) — Oee^i '?(^) nee£;\;i^'(^) • particular, the 
weight of the set A is the reliability we want to estimate, w{A) — J^aeA w{a) = 
RELs.tiG). 

Let V be the set of all paths in G from s to t. A path 7 is simply represented 
by a subset of the edges E. For every state h G H we denote by Vih) the set 
of all paths from s to t in state h. A state h € H is intact if and only if we 
have 'P{h) 0. With every path 7 S P we associate a weight ^(7) which is 
the probability that all edges on the path will be intact after the edge failure 

process, ^(7) = Uee-r 5^^) ' 

The mean number of intact paths from s to t after edge failures can thus 
be written as ^^^^-p w{'^). In Section js] we will see how this quantity can be 
efficiently calculated in acyclic graphs. 

4.2 Estimating the ratio between RELs,t{G) and the mean 
number of intact paths from s to t after edge failures 

To estimate the ratio between the reliability and the mean number of intact 
paths from s to t after edge failures we sample out of the sample space 

n = {(7, a) e 7^ X A I 7 C a}, 

where we associate the weight 111(7,0) = w{a) with every element (7,0) G 51. 
This sample space has the following interesting properties. On the one hand, 
the weight of the sample space is exactly the mean number of intact paths from 
s to t after edge failures, i.e., 

w{^)^ ^ -u;(7,a)=^ ^ u'(7, a) = ^ ^(7). 

^ y 

111(7) 

On the other hand, it is easy to sample elements of il with probability propor- 
tional to their weights by the following two-step procedure. In a first step, a 
path 7 e 7^ is drawn with probability proportional to its weight if (7). How this 
can be done efficiently will be discussed in Section [5] In the second step, all 
edges not contained in 7 are sampled corresponding to their appearance prob- 
abilities. The appeared edges of the second step together with the edges in 7 
form an intact state a G A. The tuple (7, a) is finally the sampled state. 

We will now introduce influence values for the elements in such that the 
expected influence value of a sample of is equal to w{A) /w{il). The expected 
influence value can then be estimated by a standard Monte-Carlo algorithm. 
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The approach taken in [S] was to fix for every intact state a & A an arbitrary 
s-t path 7q S 'P(a). We begin by following this approach and introduce later 
on another idea to reduce the variance of the estimator. With every sample we 
associate an infiuence value which is equal to one if the sample is of the form 
(7a, a) and equal to zero otherwise. The influence value of a sample corresponds 
therefore to the realization of a Bernoulli variable with parameter 



which is precisely the value we would like to estimate. By repeating the sampling 
procedure N times and counting the fraction of samples of the form (7a, a) we 
therefore obtain an unbiased estimator ^jv for w^Aj/w^fl) with variance 



Another unbiased estimator ijjj^ for w{A)/w(fl) with smaller variance than 
^jv can be obtained by associating an influence value 771(7,0) = 1/17^(0)1 with 
every sample (7,0), and we deflne i/jn to be simply the mean of the influence 
values of our samples. This approach is particularly interesting in our case as 
the fact of G being acyclic allows to compute the value |7^(a)| in linear time 
(see Section [s]) and thus does not increase the overall worst-case complexity of 
our algorithm. The reduction of the variance decreases the expected number 
of iterations to perform for obtaining an e — i5 approximation when applying a 
stopping criterion as presented in Section [6] 

Unfortunately, we cannot guarantee some minimal decrease of the variance 
of the estimator ip^ compared to as there are instances of reliability networks 
where we have an arbitrarily low decrease. On the other hand, it is easy to find 
instances where the variance reduces by an arbitrarily high factor. 

We finally estimate RELs^t{G) by multiplying ■(/'jv by the weight of the sam- 
pling space (i«(ri)) which can be calculated efficiently as shown in the next 
section. Algorithm [l] gives a pseudocode for an implementation of the s-t relia- 
bility estimation based on • Note that the sampling of the edges not on the 
path performed in Algorithm [T] in lines [T] to [T] can be done more efficiently as in 
general we do not need to sample all edges for calculating 1^(0)1. In Section [s] 
we will see how this part of the algorithm can be improved and show how to 
perform the remaining unspecified parts of our algorithm. More precisely, we 
will discuss the following operations, where k represents the time needed for 
generating a random number uniformly in [0, 1]: 

• Sampling in 0{mK) time a path according to line [l] of Algorithm [l] 

• Determining in 0{m) time the weight w{^) of the sampling space f2 used 
in line [T] of Algorithm [l] 

• Sampling edges not being on the initial path 7 and calculating |7^(a)| as 
needed in lines [T] to [T] of Algorithm [T] 



w({(7a,Q) en\aeA}) 
w{n) 



wiA) 
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Algorithm 1: Estimation of RELs^t{G) based on ip 



N 



Input : An acyclic reliability network G = {V,E,p) with two special 

terminals s,t G V and the number N of iterations to perform 
Output: An estimate for RELs t{G) 
a; = 0; 

for i = \ to N do 



1 

2 
3 
4 

5 
6 
7 

8 
9 
10 
11 
12 



Draw a random path 7 out of the set V with probability proportional 
to w{'y); 
a = a U 7; 

foreach e £ E\-f do 

Take a sample Xe of a Bernoulli variable with parameter q(e) to 
determine whether the edge e appears; 
if = 1 (the edge appears) then 

I a = a U e; 
end 
end 

Determine |7^(a)|; 
1 

X — X 



\V{a)\' 

14 end 

15 ipN = f ; 

16 Determine w(ri); 

17 return tpN ■ w{il.) 



5 Algorithmic details 
5.1 Sampling s-t paths 

We begin by studying how paths can be efficiently sampled according to line[l]of 
Algorithm [1] The idea is to start at the terminal s and then to construct a path 
to t by successively adding new edges. The choice of a new edge augmenting the 
current partial path is done in the following way. With every edge {v,u) G E 
we associate a weight w{v, u) which is the sum of the weights of all paths from 
V to t using edge {v,u), i.e., 

^("'")= E n^ce)- (1) 

7: path from v to t e£'y 
with (v.u)^^ 

During the path sampling method, after a partial path from s to some ver- 
tex V is constructed, we choose an outgoing edge of vertex v with probability 
proportional to the weights w. It is easy to verify that this procedure effectively 
samples a path ^ G V with probability proportional to ^(7) as desired. Further- 
more, the edge weights w can be easily computed by the following procedure. 
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We suppose without loss of generality that every edge lies on at least one 
path from s to t in G. All edges not satisfying this condition can be eliminated 
in 0{m) time in a preprocessing step. We then determine a topological order of 
the vertices. By the condition mentioned above it is clear that t will be the last 
vertex in the topological order. We go through the vertices in reverse topological 
order and determine at each step the weights w of the edges entering the current 
vertex. At the first step we look at every edge e incident to t and determine its 
weight w{e) which is equal to q{e). The weights of the other edges can then be 
determined in linear time by using the following recursive formula which follows 
from the definition of the weights w in ([T]), 

w{v,u) ^ q{v,u) ■ ^ ^(e), (2) 

eeS+(u) 

where for v €V , 5^{v) (resp. 5~{v)) denotes the set of all edges going out of v 
(resp. all edges entering v). 

5.2 Determining it>(f2) 

Another subproblcm in Algorithm [l] is the calculation of w(ri), the weight of the 
sample space. Having already calculated the weight function w, this problem 
can easily be solved by expressing w{^) in terms of w as follows. 

^ w{e) 

eeS+{s) 
—'w{e) 

5.3 Sampling edges outside the initial path and calculat- 
ing \V{a)\ 

As the sampling of the remaining edges in lines [l] to [T] of Algorithm [l] is used 
only for the calculation of |'P(a)|, we do not have to know all intact edges but 
only those on a path from s to t. One way of improving the procedure is to 
sample only edges that can be reached from s. This can be done by keeping 
track of the set of nodes C that are currently reachable from s. At the beginning, 
these are the nodes on the initial path 7. Then all edges going out of C will be 
sampled, and for every appeared edge we add its endpoint to C. This procedure 
is repeated until there are no edges outgoing from C that have not already been 
sampled. 

Finally, the determination of |'P(a)| in line [l] of Algorithm [T| can be easily 
done by an analogue technique to the one used to calculate w{^) where this 
time we work on the subgraph of G over the edges in a where we give a weight 
of one to every edge. 



eS7 
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6 Number of samples to draw 



In this section we analyze how many samples have to be drawn in our algorithm 
to obtain an e — (5 approximation. 

6.1 A priori bounds 

Using Theorem [l] we can derive the following bound. 
Theorem 3. // the number of samples N satisfies 



N>A{e- 2) In 



2\ w{n) 



then w{fl) ■ and w{fl) ■ ipN OL^e both e — S approximations of RELs.t{G) . 



In Subsection |6.3[ we discuss upper bounds for the ratio w{Q.)/w{A) which 
allow to apply Theorem [3] in practice and to formulate conditions a network has 
to satisfy under which Algorithm [T| is an FPRAS for estimating RELg t{G). 



6.2 Stopping criteria 

In practice, the use of a stopping criterion typically allows to reduce the number 
of samples to take as we can profit from information gained during the execution 
of the algorithm. Furthermore, we do not suff'er from a possibly weak upper 
bound for w{A) /w{Vl). In [2], Dagum et al. present two stopping criteria, 
applicable to our algorithm and ensuring that the result of the algorithm is an 
e — (5 approximation of REL^ ((G). In the computational results, the second 
algorithm (named Approximation Algorithm AA) is used. 



6.3 Bounding the ratio w{Q.) / w{A) 

In this subsection, we discuss upper bounds on the ratio w(f2)/u'(A). Together 
with TheoremjS] these bounds allow us to bound the number of iterations needed 
for our algorithm to deliver an e~5 approximation. The bounds discussed in this 
section are correct not only for the case of directed acyclic reliability networks 
but also for the more general case of arbitrary directed, undirected or even 
mixed reliability networks. 



Previous results Karp an Luby [Hj gave the following upper bound on the 
ratio w{^l) /w{A) (in a slightly different context). 

Theorem 4. 
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When combining the above theorem with Theorem [3] we obtain that Algo- 
rithmlllis an FPRAS if HeeBCl + ^Cs)) bounded by a polynomial in the input 
size ofthe reliability network G. In the special case of uniform edge failure 
probabilities, i.e., p(e) =p = 1 —q\/e ^ E, Theorem |4] reduces to 

and implies that if 

q = 0{log{m)/m) (4) 
then Algorithm [l] is an FPRAS for estimating RELs^t{G). 

Improved bound in the case of uniform edge failure probability We 

now give a new bound on 'w{fl)/w{A) for the case of uniform edge failure prob- 
abilities, which is sharper than ([3|, especially in the case when the reliability 
network G is not too dense and does not contain long paths from s to t. 

To quantify the sparsity of a graph we introduce the notion of edge-vertex 
bound. We say that a graph G — {V,E) has an edge- vertex bound of fi if for 
any subset of the vertices U C V we have that there are at most edges in 
the subgraph of G induced by the vertices U. The best edge-vertex bound of 
a graph can be determined in polynomial time by reduction to a flow problem 
[3]. Furthermore, Amax/'^, where Amax is the maximum degree of the vertices 
in G, can be used as a simple valid edge-vertex bound. Our new bound is given 
by the following theorem (see Appendix for a proof). 

Theorem 5. Let G — (V, E,p) be a reliability network with uniform edge failure 
probability p — I — q, edge-vertex bound fi and let Imin respectively Imax be the 
minimal and maximal length of any s-t path in G. Then we have 

^(^) < min I ^ ^ 2^+"[(«™+i"(2))(«m+i„(2)+i„„,)] ] 

w{A)- \\2-q) ' /• 

The first term of the minimum in our bound is a slight improvement over 
the bound given by ([s]). This can be seen by observing that 2/(2 ~ a;) < (1 -I- 
x) \lx S [0,1]. Contrary to bound ([S]), the first term of our bound may be 
sharp. Furthermore it is independent of the graph topology and can easily be 
generalized to non-uniform failure probabilities. 

The second term of the minimum in Theorem [5] tries to exploit some struc- 
ture of the underlying network and is particularly interesting for graphs with a 
low edge- vertex bound /i and without long paths. For example, when working 
with networks where /i is bounded by a constant and Imax — 0{^Jm log(r7i)). 
Theorem [5] implies that Algorithm [l] is an FPRAS when q = 0{y^log{m) /m), 
which is a much weaker condition than the bound given by Q . 
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7 Computational results 



In order to test our algorithm we used two random generators for creating 
directed acyclic graphs with low reliability. These generators are introduced in 
the first part of this section. In a second part, we analyze the running time 
of our algorithm on networks created by these generators with different sizes 
and reliabilities. Furthermore, the proposed algorithm is compared to a direct 
Monte-Carlo simulation. 

7.1 Test instances 

Delaunay graphs (DEL) Our generator for directed acyclic Delaunay graphs 
takes two parameters, the number of vertices n and a uniform edge intactness 
probability q. We begin by choosing n points uniformly at random in the unit 
square and consider the undirected graph given by a Delaunay triangulation of 
these points. The two terminals s and t are chosen as two vertices with maximal 
Euclidian distance. We give a linear orientation to the edges corresponding to 
the vector from s to t, i.e., an (undirected) edge {v,w} is oriented as {v,w) 
if the vector from v to w and the one from s to t have a non-negative scalar 
product, otherwise we take the orientation {w, v). Finally, all edges get uniform 
intactness probability equal to q. One can easily observe that this construction 
guarantees that every vertex lies on a path from s to t. 

Topological construction (TC) A second generator we use has three pa- 
rameters, the number of vertices n, a density parameter A G [0, 1] allowing to 
control the expected number of edges in the graph and a parameter a influenc- 
ing the intactness probabilities. We begin with an empty graph over n vertices 
V = {vi,V2, ...,«„} where Vi = s and Vn = t. The graph will be constructed 
such that {vi,V2, - ■ ■ , Vn) is a topological order of the vertices. In a first step, 
all edges of the graph are introduced, then intactness probabilities are assigned 
to the edges. 

For iG{l,2,...,n— l}we introduce an edge from Vi to i^i+i. This ensures 
that all vertices are on a path from s to t. All other possible edges will be 
present with probability A, i.e., for every i,j G {1,2, ... ,n} with i -|- 2 < j we 
add an edge {vi,Vj) with probability A. 

Finally, the intactness probability of an edge {vi,Vj) is a number chosen 
uniformly at random in the interval [0, (jUTyr^l- By choosing a < 1, edges con- 
necting topological near vertices have in general higher intactness probabilities 
than edges connecting vertices being far away from each other in the topological 
order. Therefore, smaller values for a result in less reliable networks. The value 
of a will typically be chosen in [0, 1]. 

Parameter choice The parameters were fixed in such a way that networks of 
different reliabilities were obtained for graphs with 10'', 10*, 10^ and 10^ vertices. 
We generated DEL instances for every q e {0.01,0.02, . . . , 1} for every graph 
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size n S {10^, 10^, 10^, 10^}. For the creation of TC instances the parameter A 
was always chosen such that the expected degree of every vertex is equal to ten. 
This ensures that all graphs generated with the TC generator having the same 
number of vertices also have about the same number of edges and simplifies 
the comparison of running times. For every graph size n € {10^, 10'', 10^, 10^} 
instances were generated for a £ {0.01, 0.02, . . . , 1}. 

The computational results presented in this section have been obtained on 
workstations equipped with an AMD processor 3200+ and 1GB of RAM. 



7.2 Results and interpretations 

As a first observation, we noticed that the dependence of the running time on 
e and S is essentially proportional to log(|)/e^ as predicted. We therefore fixed 
e = 0.1, 5 = 0.001 for all results presented in this section. 

Figure [l] shows the running time of the proposed algorithm as a function 
of the estimated reliability. As expected, the running time grows when larger 
reliabilities have to be estimated, as an intact state contains often several paths 
from s to t. TC instances with about the same reliability as DEL instances 
are much easier to tackle. This comes from the fact that DEL instances are 
rather locally connected, whereas most edges in TC instances were randomly 
chosen. Local connectedness has the effect that when having an intact path P 
from s to t, there is a good chance that several small subpaths of P can be 
replaced by other intact subpaths with the same start and endpoint. Every 
subpath which can be replaced in this way raises the number of intact paths 
from s to t. Furthermore, as the replaceable subpaths are small it is likely that 
large groups of them are disjoint, implying that various combinations of these 
subpath replacements yield new intact paths from s to t. Figure [l] shows that 
even instances with 10^ vertices could be solved in reasonable time as long as 
the estimated reliabihty was not too large. 
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Figure 1: Running times of the proposed algorithm in function of the 
estimated reliabilities for DEL instances (on the left) and TC instances 
(on the right). 
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For the comparison of the proposed Monte-Carlo algorithm with a direct 
Monte-Carlo approach, instances with 1000 vertices were used as most of these 
instances could have been solved in reasonable time by both algorithms. As 
well as the proposed Monte-Carlo algorithm, the direct Monte-Carlo algorithm 
was implemented by using the sampling technique explained in 5.3| allowing to 
reduce the time needed per iteration in most instances. Figure [2 shows the run- 
ning times of both algorithms on DEL and TC instances with 1000 vertices. As 
expected, the direct Monte-Carlo approach has an approximately linear depen- 
dence on the reciprocal of the estimated reliability. Figure [2] shows the strength 
of the proposed algorithm when low reliabilities have to be estimated. 
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Figure 2: Running times of the proposed algorithm (prop, algo.) and 
the direct Monte-Carlo (MC) approach for DEL instances (on the left) 
and TC instances (on the right) with 1000 vertices. 



As the running time of the direct Monte-Carlo approach depends nearly 
hnear in l/RELs,t = l/w{A) (see Figure [2]) and we expect that the proposed 
algorithm has a running time approximately proportional to wiyi) /w^A) (i.e. 
the reciprocal of the value to estimate), we expect the ratio between the run- 
ning time of the proposed algorithm and the direct Monte-Carlo approach is 
approximately linear in w(p.). This is confirmed by Figure |3] showing the ratio 
of the running time of both algorithms in function of w{^) for DEL and TC in- 
stances with 1000 vertices. It is not surprising that both algorithms need about 
the same running time when w{^) is near to one. This observation allows to 
perform a simple a priori test for deciding which algorithm is better suited for a 
particular instance. Given an instance, we first calculate w{^l) (in linear time). 
If w(ri) ^ 1 it is likely that the proposed algorithm will be faster than the direct 
Monte-Carlo approach. When w{^) 3> 1, the direct Monte-Carlo approach is 
likely to be the more efficient algorithm. 
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Figure 3: Comparison of the running time of the proposed algorithm 
and the direct Monte-Carlo approach in function of K;(ri) for networks 
with 1000 vertices. 



8 Conclusions 

An adapted version of the Monte-Carlo algorithm given by Karp and Luby in 
[9] was presented and analyzed. The new algorithm is specialized for directed 
acyclic graphs and is suited for the estimation of small reliabilities. Computa- 
tional results show the successful application of the proposed algorithm on two 
types of randomly generated large-scale instances and its advantage compared 
to the direct Monte-Carlo approach when very small reliabilities have to be es- 
timated. Previous algorithms for accurate estimation of s-t reliability were only 
applicable on either very small instances or on a very restricted class of initial 
networks. For the case of uniform edge failure probabilities, a worst-case bound 
on the number of samples to be drawn was given that sharpens a bound pre- 
sented in and is significantly stronger in the case of relatively sparse graphs 
without long paths from s to t. 

One important open question in this domain is if there exists an FPRAS for 
estimating s-t reliability in directed acyclic graphs. It would be interesting to 
find algorithms allowing to tackle instances efficiently that cannot be solved in 
reasonable time by our algorithm or the direct Monte-Carlo approach. Another 
point is the generalization of the upper bound for w{^) lw{A) for the case of 
non-uniform failure probabilities. Additionally for the case of general (not nec- 
essarily acyclic) networks there seem to be no practically efficient algorithms at 
the moment for the estimation of low reliabilities on large instances. 
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Appendix 

Proof of Theorem [5] 

Let (7, a) e be a sample drawn according to lines [T][T] of Algorithm [T| We 
define the influence value m(7, a) of the sample as in Section |4] as a random 
variable ^ 

m-.n — >R, m(7, a) -y . 

As discussed in Section |4] the influence value of a sample from il can be 
used as an unbiased estimator for w{A)/w{n) as we have 

E[™(7,a)] = ^. 

By the above equation, the reciprocal of a lower bound on E[m] is an upper 
bound on w{fl)/w{A). We will therefore deduce the bound given in Theorem |5] 
by deriving a lower bound bound on w{A)/w{fl). 

Let 0.1 k be the subset of all elements of the sample space where the initial 
chosen path has length I and the total number of appeared edges is I + k, i.e., 

ni^k^{ii,a) en\ I7I = ;,|a| = ; + fc} . 

It is clear that all elements u g fli ^ appear with equal probability and can 
be sampled by the method described in Algorithm [2] This method is introduced 
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just for theoretical analysis. 



Algorithm 2: Sampling uniformly from Q 



Lk 



1 Choose uniformly a path 7 G T' with length I; 

2 a = 7; 

3 R = E\r, 

4 for i = 1 to fc do 

5 Choose an edge e i? uniformly at random; 

6 a — a U Bi] 

7 R= R\ef, 

8 end 

9 return (7, a) 



With every edge added in the for-loop of Algorithm |2] we associate a 
multiplicity Afg. G {1/2,1} equal to 1/2 if both endpoints of are saturated 
by edges in 7 U {e^ | 1 < j < « — 1} and equal to 1 otherwise. Intuitively, the 
multiplicity is a measure for the influence of an added edge in Algorithm [2] on 
the ratio w{A)/w(fl). The following lemma formalizes this intuition. 

Lemma 1. Let (7, a) e f2;^fc be an intact state constructed as described in 
Algorithm^ Then we have 

Before proving this lemma, we discuss the link between the lemma and 
Algorithm [2j In general there are multiple ways to get an intact state a G A 
with Algorithm[2] Depending on how the state a was obtained, the multiplicities 
associated with the edges in a are different. As Lemma[l]is true for every possible 
way of obtaining state a by Algorithm |2] it is applicable even to intact states 
that were not constructed through Algorithm |2] One just has to fix a possible 
way how the state a could have been constructed by Algorithm [2] i.e., an s-t 
path 7 C a has to be fixed as well as an order for the edges in a \ 7 specifying in 
which sequence those edges were chosen in Algorithm [2j The multiplicities can 
then be calculated with respect to this order, and Lemma [T| can be applied. 

Proof of Lemma^ Let (7, a) G fli^k be an intact state constructed as described 
in Algorithm|2] We define F — {ei, 62, . . . , e^} = a \ 7 to be the set of all edges 
added to the initial path 7 during the construction of a. Let F'i^F^ be the 
following partitioning of the edges in a. 

F"^={eeF\M, = ^} 

= 7 U {e e F I A//e = 1} 

We prove the following statement, which immediately implies Lemma [l] 

Proposition. For every set H C F^ there exists at most one s-t path in the 
state a that contains all edges of H and none of F^ \ H. 
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Lemma [T] follows from the above proposition by the following observation. 
The proposition implies that there are at most as many different s-t paths in a 
as there are subsets of F2 . We therefore have 



which finally implies 



\ria)\<2 



F2 



\'Pia)\ 



> 



\F2 



It remains to prove the proposition. Let H C and suppose that we have 
two different s-t paths 71 , 72 ^ a such that both contain all edges of H and none 
of F2 \ H. This implies that their symmetric difference 71A72 C a contains a 
cycle consisting only of edges in F^, thus contradicting the fact that F^ does 
not contain cycles (any cycle in a contains at least one element of Fs). 

□ 

With the aid of Lemma [T] we prove the following intermediate result. 

Lemma 2. Let I £ N,k E {0} U N with I + k < m and uj = (7, a) be a random 
sample from the sample space according to Algorithm^ We have 



E 



\na)\ 



(7, a) e flik 



-k{k+l) 



Proof of Lemma^ Let uj' — (7', a') be a sample corresponding to a result of 
Algorithm [2] for the given I and k. Conditioned on ^li^k, w has therefore the 

, Afgfc be the random variables corre- 



same distribution as 



Let Me^ , Me^ , . 



sponding to the multiplicities of the edges in a' \ 7'. By Lemma [T] we have 



E 



\V{a) 



(7, a) g fli^k 



= E 



\na')\ 



> E 



(5) 



Let i € {1, . . . , fc}. Observe that independently of which path 7' and which 
edges ei, . . . , efc were chosen, we have that at most ^ + 1 + 2(i — 1) vertices in 
G are saturated by the edges in 7' U {ei, . . . , Ci-i} {I + 1 vertices are saturated 
through 7' and every additional edge saturates at most two new vertices in G). 
Let Vi be the vertices in G which are saturated by 7' U {ei, . . . , e^-i}. 

In Algorithm |2] the edge is chosen uniformly at random from the remain- 
ing edges ii^\(7'U{ei, . . . , Cfc}). As G is sparse with edge- vertex bound fi we have 
that at most n\V'^\ < fi{l + l + 2{i — l)) edges have both endpoints in Vi. Further- 
more, « — 1 of these edges were already chosen. We therefore have the following 
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stochastic inequality (which is true for any reaHzation of 7', Mg^, . . . , Afe. ^): 
P fM. = i I 7', Afe. . . . , Afe._.) < min (ir'^^l - - 1) 



< mm < 1 , 

[ m 

r n(2i + l-l) 

< iJimi 1, ^ — '- 

[ m 

The stochastic inequahty above allows to give a simple bound on the follow- 
ing conditional expectation: 

E[Me, I y, Me,, ... , Me._J = ^ P (^A^, = ^ | 7', , . . . , Me,_,^ + 

P(Afe. = 1 |7',Me,,...,Me._J 
= l-^p(^Me, =.^|7',Me,,...,Afe._, 

> max < - , 1 



.2' 2m 

Applying the above inequality, the expectation in ^ can be developed as 



E 



>nmax|-,l 1 

>nmaxJi,Q) 



fi(2i + i-l) 



The last inequality comes from the fact that (1/4)^ < 1 — a; for a; g [0,1/2] 
and (1/4)^ < 1/2 for x > 1/2. Developing further, we finally get the result of 
Lemma [2l 



riniax \ J, f\) I > max i 2-^2".S 



i=l 



H(2i + i-l) 



>max{2-^2-^'=('=+')} 



□ 



Beginning with the result of Lemma |2] we now prove Theorem |5] by first 
weakening and then eliminating the conditioning on ^i^k- Let fli be the set of 
all elements of the sample space where the initial chosen path has length I, 
i.e.. 
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Let K be the random variable corresponding to the number of edges that 
appeared additionally to the ones of the initial path, when drawing an element 
out of fi; . Note that K is binomially distributed as 



Using Lemma [2] we get 



E 



v\na)\ 



(7, a) e 



K ^ Bin(m — l,q) 



> E 



max {2-^, 2- ^^(^'+')} 
>max{E[2-^] , E [2-^^^^+')] } 



1 - 



m—l 



,E 



(6) 



By replacing I by Imin in the first term of the above maximum we get the 
first part of the inequality in Theorem [5] i.e., 



wiA) 



E 



l^(«)l 



< 1- 



-(m-i„i„) 



2-q 



The remaining part of Theorem [s] will be shown by developing the second 
term of ^ further. We will use the inequality 

E[2-^^(^+')] > 2-^°("+')p[i^ < a] , 

which is true for any a S E and a result shown by Hamza [5] stating that the 
distance between the median and the mean of a binomial random variable is at 
most ln(2). Therefore, when choosing a = qrn + ln{2), we have P[K < a] > 1/2 
and get the following result: 

gj2-™-'^(-'^+')] > 2"l~TS:(«('"~')+'"(2))(9(m-0+In(2)+0 
> 2-l-J^(9™+ln(2))(9™+ln(2)+0 

This implies the second part of the inequality in Theorem |5] as 



w{A) 



E 



^ 2l+^(9™+ln(2))(9"'+ln(2)+/n,ax) 



which completes the proof of Theorem [5j 
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